Acute Changes in CEWL as a Function of Temperature: Data Analysis

Calvin Davis, Savannah Weaver

2022

Packages

if (!require("tidyverse")) install.packages("tidyverse") 
library("tidyverse") #tidyr, ggplot, dplyr, etc
if (!require("ggpubr")) install.packages("ggpubr") 
library("ggpubr")
if (!require("lme4")) install.packages("lme4") 
library("lme4") # mixed-effect models
if (!require("lmerTest")) install.packages("lmerTest") 
library("lmerTest") # p-values
if (!require("car")) install.packages("car") 
library("car") # test for VIF
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans')
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # model export
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown format

Background

This code walks through the statistical models and figures for a study on the rate of change of cutaneous evaporative water loss (CEWL) in response to temperature change in Western Fence Lizards (Sceloporus occidentalis). Data was collected and analyzed by Calvin Davis and Savannah Weaver, and the project was advised by Dr. Emily Taylor at California Polytechnic State University. The sections with an asterisk are the stats and figures presented in the publication.

Load Data

Get the RDS files created in the data wrangling Rmd:

lizard_tmt_info <- read_rds("Data/collection_dat_formatted.RDS")
temp_time_series <- read_rds("Data/temp_time_series.RDS") %>% 
  dplyr::filter(lizard_ID != 410)
CEWL_time_series <- read_rds("Data/CEWL_time_series.RDS") # not needed???
unique(CEWL_time_series$date)
## [1] "2021-10-05" "2021-10-11" "2021-10-12" "2021-10-18" "2021-10-19"
by_minute_indiv <- temp_time_series %>% 
  mutate(time_min_block = floor(time_min)) %>% 
  group_by(date, lizard_ID, treatment, time_min_block) %>% 
  summarise(CEWL_per_min = mean(CEWL_g_m2h),
            clo_temp_per_min = mean(calibrated_cloacal_temp),
            dors_temp_per_min = mean(calibrated_dorsal_temp)) 
## `summarise()` has grouped output by 'date', 'lizard_ID', 'treatment'. You can
## override using the `.groups` argument.
by_minute_avg <- by_minute_indiv %>% 
  group_by(treatment, time_min_block) %>% 
  summarise(CEWL_min_avg = mean(CEWL_per_min),
            clo_temp_min_avg = mean(clo_temp_per_min),
            dors_temp_min_avg = mean(dors_temp_per_min))
## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.
by_temp_minute_avg <- by_minute_indiv %>% 
  mutate(clo_temp_chunk = floor(clo_temp_per_min)) %>% 
  group_by(treatment, clo_temp_chunk) %>% 
  summarise(CEWL_min_avg = mean(CEWL_per_min))
## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.

Stats & Models

* Temp Correlation

Check the correlation between dorsal and cloacal temps:

temp_corr <- lm(data = by_minute_avg, 
                dors_temp_min_avg ~ clo_temp_min_avg)
summary(temp_corr)
## 
## Call:
## lm(formula = dors_temp_min_avg ~ clo_temp_min_avg, data = by_minute_avg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3083 -0.5866  0.1572  0.2957  1.6472 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.69880    0.65881   10.17 5.19e-13 ***
## clo_temp_min_avg  0.75694    0.02499   30.29  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7538 on 43 degrees of freedom
## Multiple R-squared:  0.9552, Adjusted R-squared:  0.9542 
## F-statistic: 917.7 on 1 and 43 DF,  p-value: < 2.2e-16

CEWL (raw) ~ temp LMM

First, test linear effects of temperature:

CEWL_LMM_02a <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             clo_temp_per_min *
                             treatment +
                             (1|date/lizard_ID))

CEWL_LMM_02b <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             dors_temp_per_min *
                             treatment +
                             (1|date/lizard_ID))

Compare to the inclusion of polynomial factors:

CEWL_LMM_03a <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             poly(clo_temp_per_min, 2)*treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_03a, CEWL_LMM_02a)
## refitting model(s) with ML (instead of REML)
## Data: by_minute_indiv
## Models:
## CEWL_LMM_02a: CEWL_per_min ~ clo_temp_per_min * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_per_min ~ poly(clo_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_02a    9 1786.0 1823.1 -883.98   1768.0                         
## CEWL_LMM_03a   12 1690.4 1740.0 -833.20   1666.4 101.56  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CEWL_LMM_03b <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             poly(dors_temp_per_min, 2)*treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_03b, CEWL_LMM_02b)
## refitting model(s) with ML (instead of REML)
## Data: by_minute_indiv
## Models:
## CEWL_LMM_02b: CEWL_per_min ~ dors_temp_per_min * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_per_min ~ poly(dors_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_02b    9 1809.6 1846.8 -895.80   1791.6                         
## CEWL_LMM_03b   12 1724.3 1773.8 -850.15   1700.3 91.304  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The poly models have significantly better explanatory power than the linear models.

Also try log transformation to make sure a polynomial effect is best. But, cannot compare a log-log model because then the response variable would be different.

CEWL_LMM_04a <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             log(clo_temp_per_min)*treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_04a, CEWL_LMM_02a)
## refitting model(s) with ML (instead of REML)
## Data: by_minute_indiv
## Models:
## CEWL_LMM_04a: CEWL_per_min ~ log(clo_temp_per_min) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_02a: CEWL_per_min ~ clo_temp_per_min * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## CEWL_LMM_04a    9 1801.4 1838.6 -891.72   1783.4                     
## CEWL_LMM_02a    9 1786.0 1823.1 -883.98   1768.0 15.475  0
anova(CEWL_LMM_04a, CEWL_LMM_03a)
## refitting model(s) with ML (instead of REML)
## Data: by_minute_indiv
## Models:
## CEWL_LMM_04a: CEWL_per_min ~ log(clo_temp_per_min) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_per_min ~ poly(clo_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_04a    9 1801.4 1838.6 -891.72   1783.4                         
## CEWL_LMM_03a   12 1690.4 1740.0 -833.20   1666.4 117.03  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CEWL_LMM_04b <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             log(dors_temp_per_min)*treatment +
                             treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_04b, CEWL_LMM_02b)
## refitting model(s) with ML (instead of REML)
## Data: by_minute_indiv
## Models:
## CEWL_LMM_04b: CEWL_per_min ~ log(dors_temp_per_min) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_02b: CEWL_per_min ~ dors_temp_per_min * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## CEWL_LMM_04b    9 1825.5 1862.6 -903.73   1807.5                     
## CEWL_LMM_02b    9 1809.6 1846.8 -895.80   1791.6 15.858  0
anova(CEWL_LMM_04b, CEWL_LMM_03b)
## refitting model(s) with ML (instead of REML)
## Data: by_minute_indiv
## Models:
## CEWL_LMM_04b: CEWL_per_min ~ log(dors_temp_per_min) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_per_min ~ poly(dors_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_04b    9 1825.5 1862.6 -903.73   1807.5                         
## CEWL_LMM_03b   12 1724.3 1773.8 -850.15   1700.3 107.16  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The logarithmic models are NOT better than the linear or polynomial ones. The poly models are best and the linear and log models are equally less-good.

Check assumptions:

plot(CEWL_LMM_03a)

plot(CEWL_LMM_03b)

vif(CEWL_LMM_03a)
##                                             GVIF Df GVIF^(1/(2*Df))
## poly(clo_temp_per_min, 2)           788782.98078  2       29.801586
## treatment                                5.82529  2        1.553565
## poly(clo_temp_per_min, 2):treatment 787225.05195  4        5.457734
vif(CEWL_LMM_03b)
##                                              GVIF Df GVIF^(1/(2*Df))
## poly(dors_temp_per_min, 2)           4.192428e+05  2       25.445817
## treatment                            5.645315e+00  2        1.541424
## poly(dors_temp_per_min, 2):treatment 4.183579e+05  4        5.043053

There’s no fanning or different pattern. VIF is negligible.

Re-run with p-values:

CEWL_LMM_besta <- lmerTest::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             treatment * poly(clo_temp_per_min, 2) +
                             (1|date/lizard_ID))
summary(CEWL_LMM_besta)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## CEWL_per_min ~ treatment * poly(clo_temp_per_min, 2) + (1 | date/lizard_ID)
##    Data: by_minute_indiv
## 
## REML criterion at convergence: 1625.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3166 -0.5957 -0.0968  0.5423  5.8100 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 12.210   3.494   
##  date           (Intercept) 24.154   4.915   
##  Residual                    1.578   1.256   
## Number of obs: 459, groups:  lizard_ID:date, 32; date, 5
## 
## Fixed effects:
##                                              Estimate Std. Error        df
## (Intercept)                                   -0.4704     3.7914   29.3064
## treatmentCooling                               8.2628     3.2867  261.5997
## treatmentHeating                              17.9752     3.2917  244.1138
## poly(clo_temp_per_min, 2)1                   139.9853    50.9856  441.7664
## poly(clo_temp_per_min, 2)2                  -502.4452    82.9294  445.3354
## treatmentCooling:poly(clo_temp_per_min, 2)1  -99.1502    51.0336  441.7445
## treatmentHeating:poly(clo_temp_per_min, 2)1 -114.1937    51.0301  441.7380
## treatmentCooling:poly(clo_temp_per_min, 2)2  519.6856    82.9658  445.3323
## treatmentHeating:poly(clo_temp_per_min, 2)2  490.2494    82.9555  445.3375
##                                             t value Pr(>|t|)    
## (Intercept)                                  -0.124  0.90211    
## treatmentCooling                              2.514  0.01254 *  
## treatmentHeating                              5.461 1.17e-07 ***
## poly(clo_temp_per_min, 2)1                    2.746  0.00629 ** 
## poly(clo_temp_per_min, 2)2                   -6.059 2.92e-09 ***
## treatmentCooling:poly(clo_temp_per_min, 2)1  -1.943  0.05267 .  
## treatmentHeating:poly(clo_temp_per_min, 2)1  -2.238  0.02573 *  
## treatmentCooling:poly(clo_temp_per_min, 2)2   6.264 8.86e-10 ***
## treatmentHeating:poly(clo_temp_per_min, 2)2   5.910 6.82e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC trtmnH p(___,2)1 p(___,2)2 tC:(___,2)1 tH:(___,2)1
## tretmntClng -0.769                                                          
## tretmntHtng -0.764  0.882                                                   
## ply(___,2)1 -0.709  0.820  0.814                                            
## ply(___,2)2  0.760 -0.882 -0.872 -0.878                                     
## tC:(___,2)1  0.708 -0.818 -0.814 -0.999     0.877                           
## tH:(___,2)1  0.708 -0.819 -0.814 -0.999     0.878     0.998                 
## tC:(___,2)2 -0.760  0.881  0.872  0.878    -1.000    -0.877      -0.877     
## tH:(___,2)2 -0.760  0.882  0.872  0.878    -1.000    -0.877      -0.878     
##             tC:(___,2)2
## tretmntClng            
## tretmntHtng            
## ply(___,2)1            
## ply(___,2)2            
## tC:(___,2)1            
## tH:(___,2)1            
## tC:(___,2)2            
## tH:(___,2)2  0.999
CEWL_clo_anova <- data.frame(anova(CEWL_LMM_besta, type = "2", ddf = "Kenward-Roger")) %>% 
  mutate(variable = row.names(.))

CEWL_LMM_bestb <- lmerTest::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             treatment * poly(dors_temp_per_min, 2) +
                             (1|date/lizard_ID))
summary(CEWL_LMM_bestb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_per_min ~ treatment * poly(dors_temp_per_min, 2) + (1 |  
##     date/lizard_ID)
##    Data: by_minute_indiv
## 
## REML criterion at convergence: 1660.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4453 -0.6089 -0.0814  0.5142  5.9945 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 10.295   3.209   
##  date           (Intercept) 22.052   4.696   
##  Residual                    1.728   1.315   
## Number of obs: 459, groups:  lizard_ID:date, 32; date, 5
## 
## Fixed effects:
##                                              Estimate Std. Error       df
## (Intercept)                                     2.194      3.495   25.854
## treatmentCooling                                4.972      2.973  262.431
## treatmentHeating                               15.455      2.988  246.451
## poly(dors_temp_per_min, 2)1                   179.875     54.756  444.012
## poly(dors_temp_per_min, 2)2                  -410.895     72.930  441.411
## treatmentCooling:poly(dors_temp_per_min, 2)1 -149.678     54.791  443.997
## treatmentHeating:poly(dors_temp_per_min, 2)1 -158.197     54.802  443.991
## treatmentCooling:poly(dors_temp_per_min, 2)2  428.187     72.966  441.398
## treatmentHeating:poly(dors_temp_per_min, 2)2  403.389     72.965  441.425
##                                              t value Pr(>|t|)    
## (Intercept)                                    0.628  0.53572    
## treatmentCooling                               1.672  0.09565 .  
## treatmentHeating                               5.173 4.78e-07 ***
## poly(dors_temp_per_min, 2)1                    3.285  0.00110 ** 
## poly(dors_temp_per_min, 2)2                   -5.634 3.14e-08 ***
## treatmentCooling:poly(dors_temp_per_min, 2)1  -2.732  0.00655 ** 
## treatmentHeating:poly(dors_temp_per_min, 2)1  -2.887  0.00408 ** 
## treatmentCooling:poly(dors_temp_per_min, 2)2   5.868 8.65e-09 ***
## treatmentHeating:poly(dors_temp_per_min, 2)2   5.529 5.54e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC trtmnH p(___,2)1 p(___,2)2 tC:(___,2)1 tH:(___,2)1
## tretmntClng -0.753                                                          
## tretmntHtng -0.748  0.878                                                   
## ply(___,2)1 -0.721  0.849  0.843                                            
## ply(___,2)2  0.743 -0.877 -0.868 -0.923                                     
## tC:(___,2)1  0.721 -0.848 -0.843 -0.999     0.923                           
## tH:(___,2)1  0.720 -0.848 -0.843 -0.999     0.923     0.999                 
## tC:(___,2)2 -0.743  0.876  0.867  0.923    -1.000    -0.922      -0.922     
## tH:(___,2)2 -0.743  0.876  0.867  0.923    -0.999    -0.922      -0.922     
##             tC:(___,2)2
## tretmntClng            
## tretmntHtng            
## ply(___,2)1            
## ply(___,2)2            
## tC:(___,2)1            
## tH:(___,2)1            
## tC:(___,2)2            
## tH:(___,2)2  0.999
CEWL_dors_anova <- data.frame(anova(CEWL_LMM_bestb, type = "2", ddf = "Kenward-Roger")) %>%
  mutate(variable = row.names(.))

# double check sample sizes
temp_time_series %>%
  group_by(lizard_ID, treatment) %>%
  summarise(n()) %>%
  group_by(treatment) %>%
  summarise(n())
## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.
## # A tibble: 3 × 2
##   treatment `n()`
##   <fct>     <int>
## 1 Control      11
## 2 Cooling      11
## 3 Heating      10

These are the best models for CEWL ~ temperature.

Export:

# model anovas
CEWL_temp_F_stats <- CEWL_clo_anova %>%
  rbind(CEWL_dors_anova) %>% 
  mutate(partial_sum_of_squares = round(Sum.Sq, 0), 
         F_only = round(F.value, 2),
         DenDF = round(DenDF, 0),
         F_statistic = paste(F_only, " (", NumDF, ",", DenDF, ")", sep = ""),
         p_value = (Pr..F.)) %>%
  dplyr::select(variable,
                partial_sum_of_squares,
                F_statistic, 
                p_value)
write.csv(CEWL_temp_F_stats, 
          "./Results_Statistics/CEWL_temp_F_stats.csv")

# model summaries
broom.mixed::tidy(CEWL_LMM_besta) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_clotemp.csv")

broom.mixed::tidy(CEWL_LMM_bestb) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_dorstemp.csv")

CEWL (raw) ~ time LMM

NOTE: Start times vary from 97-170 seconds after starting time series.

Make a polynomial model first:

rates_LMM_01 <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             poly(time_min_block, 2) * treatment +
                             (1|date/lizard_ID))
summary(rates_LMM_01)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## CEWL_per_min ~ poly(time_min_block, 2) * treatment + (1 | date/lizard_ID)
##    Data: by_minute_indiv
## 
## REML criterion at convergence: 1630.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6540 -0.5335  0.0275  0.4644  5.8089 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.089   3.015   
##  date           (Intercept) 20.204   4.495   
##  Residual                    1.581   1.258   
## Number of obs: 459, groups:  lizard_ID:date, 32; date, 5
## 
## Fixed effects:
##                                           Estimate Std. Error t value
## (Intercept)                                 15.273      2.211   6.909
## poly(time_min_block, 2)1                   -17.822      2.119  -8.409
## poly(time_min_block, 2)2                     6.503      2.145   3.032
## treatmentCooling                            -8.441      1.312  -6.435
## treatmentHeating                             2.310      1.361   1.697
## poly(time_min_block, 2)1:treatmentCooling  -17.283      3.107  -5.562
## poly(time_min_block, 2)2:treatmentCooling   15.775      3.074   5.132
## poly(time_min_block, 2)1:treatmentHeating   41.968      3.080  13.628
## poly(time_min_block, 2)2:treatmentHeating  -14.173      3.082  -4.598
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(t__,2)1 -0.001                                                        
## ply(t__,2)2  0.000  0.001                                                 
## tretmntClng -0.292  0.002     0.000                                       
## tretmntHtng -0.286  0.002     0.000     0.466                             
## pl(__,2)1:C  0.001 -0.682     0.000     0.005 -0.001                      
## pl(__,2)2:C  0.000  0.000    -0.698     0.002  0.000  0.025               
## pl(__,2)1:H  0.001 -0.688     0.000    -0.002 -0.001  0.469      0.000    
## pl(__,2)2:H  0.000  0.000    -0.696     0.000  0.000  0.000      0.486    
##             p(__,2)1:H
## ply(t__,2)1           
## ply(t__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.008

Also make a linear model version:

rates_LMM_02 <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             time_min_block * treatment +
                             (1|date/lizard_ID))

anova(rates_LMM_02, rates_LMM_01)
## refitting model(s) with ML (instead of REML)
## Data: by_minute_indiv
## Models:
## rates_LMM_02: CEWL_per_min ~ time_min_block * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_per_min ~ poly(time_min_block, 2) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## rates_LMM_02    9 1786.6 1823.8 -884.31   1768.6                         
## rates_LMM_01   12 1682.6 1732.2 -829.30   1658.6 110.02  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Polynomial is much better than linear.

Also compare to log version:

rates_LMM_03 <- lme4::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             log(time_min_block) * treatment +
                             (1|date/lizard_ID))

rates_LMM_03_loglog <- lme4::lmer(data = by_minute_indiv,
                          log(CEWL_per_min) ~ 
                             log(time_min_block) * treatment +
                             (1|date/lizard_ID))

The polynomial model is best (slightly better than log and much better than lm options).

Check linear regression assumptions:

plot(rates_LMM_01)

vif(rates_LMM_01)
##                                       GVIF Df GVIF^(1/(2*Df))
## poly(time_min_block, 2)           8.011443  2        1.682394
## treatment                         1.000098  2        1.000024
## poly(time_min_block, 2):treatment 8.011931  4        1.297081

re-run with p-values:

rates_LMM_best <- lmerTest::lmer(data = by_minute_indiv,
                          CEWL_per_min ~ 
                             poly(time_min_block, 2) * treatment +
                             (1|date/lizard_ID))
summary(rates_LMM_best)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## CEWL_per_min ~ poly(time_min_block, 2) * treatment + (1 | date/lizard_ID)
##    Data: by_minute_indiv
## 
## REML criterion at convergence: 1630.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6540 -0.5335  0.0275  0.4644  5.8089 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.089   3.015   
##  date           (Intercept) 20.204   4.495   
##  Residual                    1.581   1.258   
## Number of obs: 459, groups:  lizard_ID:date, 32; date, 5
## 
## Fixed effects:
##                                           Estimate Std. Error      df t value
## (Intercept)                                 15.273      2.211   5.075   6.909
## poly(time_min_block, 2)1                   -17.822      2.119 421.033  -8.409
## poly(time_min_block, 2)2                     6.503      2.145 421.043   3.032
## treatmentCooling                            -8.441      1.312  25.139  -6.435
## treatmentHeating                             2.310      1.361  25.208   1.697
## poly(time_min_block, 2)1:treatmentCooling  -17.283      3.107 421.456  -5.562
## poly(time_min_block, 2)2:treatmentCooling   15.775      3.074 421.092   5.132
## poly(time_min_block, 2)1:treatmentHeating   41.968      3.080 421.102  13.628
## poly(time_min_block, 2)2:treatmentHeating  -14.173      3.082 421.059  -4.598
##                                           Pr(>|t|)    
## (Intercept)                               0.000916 ***
## poly(time_min_block, 2)1                  6.47e-16 ***
## poly(time_min_block, 2)2                  0.002582 ** 
## treatmentCooling                          9.49e-07 ***
## treatmentHeating                          0.101948    
## poly(time_min_block, 2)1:treatmentCooling 4.74e-08 ***
## poly(time_min_block, 2)2:treatmentCooling 4.38e-07 ***
## poly(time_min_block, 2)1:treatmentHeating  < 2e-16 ***
## poly(time_min_block, 2)2:treatmentHeating 5.64e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(t__,2)1 -0.001                                                        
## ply(t__,2)2  0.000  0.001                                                 
## tretmntClng -0.292  0.002     0.000                                       
## tretmntHtng -0.286  0.002     0.000     0.466                             
## pl(__,2)1:C  0.001 -0.682     0.000     0.005 -0.001                      
## pl(__,2)2:C  0.000  0.000    -0.698     0.002  0.000  0.025               
## pl(__,2)1:H  0.001 -0.688     0.000    -0.002 -0.001  0.469      0.000    
## pl(__,2)2:H  0.000  0.000    -0.696     0.000  0.000  0.000      0.486    
##             p(__,2)1:H
## ply(t__,2)1           
## ply(t__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.008
anova(rates_LMM_best, type = "2", ddf = "Kenward-Roger")
## Type II Analysis of Variance Table with Kenward-Roger's method
##                                   Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
## poly(time_min_block, 2)           160.58  80.289     2 421.20  50.768 < 2.2e-16
## treatment                         108.67  54.337     2  25.21  34.358 6.314e-08
## poly(time_min_block, 2):treatment 736.17 184.044     4 421.19 116.374 < 2.2e-16
##                                      
## poly(time_min_block, 2)           ***
## treatment                         ***
## poly(time_min_block, 2):treatment ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# double check sample sizes
CEWL_time_series %>%
  group_by(lizard_ID, treatment) %>%
  summarise(n()) %>%
  group_by(treatment) %>%
  summarise(n())
## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.
## # A tibble: 3 × 2
##   treatment `n()`
##   <fct>     <int>
## 1 Control      12
## 2 Cooling      12
## 3 Heating      13

Export:

broom.mixed::tidy(rates_LMM_best) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_time.csv")

Rate of Change

* Tb ~ Time

lizard_temp_time_lm_log <- lm(data = temp_time_series, 
                          log(calibrated_cloacal_temp) ~ log(time_min)*lizard_ID)
plot(lizard_temp_time_lm_log)

Now get the trend/slope for each individual and make into a variable in df.

lizard_coeffs_temp <- data.frame(emtrends(lizard_temp_time_lm_log, 
                                          "lizard_ID", var = "time_min")) %>%
   left_join(lizard_tmt_info, 
             by = c("lizard_ID" = "individual_ID")) %>%
   dplyr::select(lizard_ID, rate_change = time_min.trend, 
                 SE, df, lower.CL, upper.CL,
                 date,
                 treatment
                 ) %>%
   dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs_temp)
##    lizard_ID   rate_change               SE                  df       
##  401    : 1   Min.   :-0.0545701   Min.   :0.0001623   Min.   :24381  
##  402    : 1   1st Qu.:-0.0324598   1st Qu.:0.0001752   1st Qu.:24381  
##  404    : 1   Median : 0.0021493   Median :0.0001779   Median :24381  
##  406    : 1   Mean   :-0.0001463   Mean   :0.0001863   Mean   :24381  
##  407    : 1   3rd Qu.: 0.0346028   3rd Qu.:0.0001876   3rd Qu.:24381  
##  408    : 1   Max.   : 0.0440239   Max.   :0.0003102   Max.   :24381  
##  (Other):26                                                           
##     lower.CL             upper.CL               date              treatment 
##  Min.   :-0.0549235   Min.   :-0.0542167   Min.   :2021-10-05   Control:11  
##  1st Qu.:-0.0328089   1st Qu.:-0.0321107   1st Qu.:2021-10-11   Cooling:11  
##  Median : 0.0018033   Median : 0.0024954   Median :2021-10-12   Heating:10  
##  Mean   :-0.0005115   Mean   : 0.0002189   Mean   :2021-10-13               
##  3rd Qu.: 0.0342589   3rd Qu.: 0.0349467   3rd Qu.:2021-10-18               
##  Max.   : 0.0436800   Max.   : 0.0443679   Max.   :2021-10-19               
##                                                                             
##      date_factor
##  2021-10-05:6   
##  2021-10-11:6   
##  2021-10-12:7   
##  2021-10-18:7   
##  2021-10-19:6   
##                 
## 
length(unique(lizard_coeffs_temp$lizard_ID))
## [1] 32

then get average rate of change by treatment group:

avg_temp_change <- lm(data = lizard_coeffs_temp, rate_change ~ treatment)
summary(avg_temp_change)
## 
## Call:
## lm(formula = rate_change ~ treatment, data = lizard_coeffs_temp)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0170679 -0.0020104  0.0005646  0.0021987  0.0134091 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.002336   0.001653   1.413    0.168    
## treatmentCooling -0.039838   0.002338 -17.042  < 2e-16 ***
## treatmentHeating  0.035878   0.002395  14.978 3.49e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005482 on 29 degrees of freedom
## Multiple R-squared:  0.9719, Adjusted R-squared:   0.97 
## F-statistic: 501.3 on 2 and 29 DF,  p-value: < 2.2e-16
avg_temp_change_coeff <- data.frame(emmeans(avg_temp_change, "treatment"))
avg_temp_change_coeff
##   treatment       emmean          SE df    lower.CL     upper.CL
## 1   Control  0.002336053 0.001652941 29 -0.00104459  0.005716696
## 2   Cooling -0.037502169 0.001652941 29 -0.04088281 -0.034121526
## 3   Heating  0.038214475 0.001733619 29  0.03466883  0.041760123
avg_temp_change_pairs <- data.frame(pairs(emmeans(avg_temp_change, "treatment")))
avg_temp_change_pairs
##            contrast    estimate          SE df   t.ratio      p.value
## 1 Control - Cooling  0.03983822 0.002337611 29  17.04228 3.885781e-15
## 2 Control - Heating -0.03587842 0.002395338 29 -14.97844 1.454392e-14
## 3 Cooling - Heating -0.07571664 0.002395338 29 -31.61000 3.552714e-15

and finally, t-tests of rate of change per tmt group:

lizard_coeffs_temp_control <- lizard_coeffs_temp %>% 
  dplyr::filter(treatment == "Control")
t.test(lizard_coeffs_temp_control$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_temp_control$rate_change
## t = 6.5738, df = 10, p-value = 6.282e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.001544266 0.003127840
## sample estimates:
##   mean of x 
## 0.002336053
lizard_coeffs_temp_cool <- lizard_coeffs_temp %>% 
  dplyr::filter(treatment == "Cooling")
t.test(lizard_coeffs_temp_cool$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_temp_cool$rate_change
## t = -14.558, df = 10, p-value = 4.66e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.04324210 -0.03176224
## sample estimates:
##   mean of x 
## -0.03750217
lizard_coeffs_temp_heat <- lizard_coeffs_temp %>% 
  dplyr::filter(treatment == "Heating")
t.test(lizard_coeffs_temp_heat$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_temp_heat$rate_change
## t = 32.083, df = 9, p-value = 1.364e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.03552000 0.04090895
## sample estimates:
##  mean of x 
## 0.03821447

* CEWL ~ Time

We want to know why each lizard had a different rate of change of CEWL during the 15-minute temperature treatment and measurement period.

First, calculate rate of change for each lizard:

lizard_CEWL_time_lm_log <- lm(data = temp_time_series,
                          log(CEWL_g_m2h) ~ 
                             log(time_min)*lizard_ID)
plot(lizard_CEWL_time_lm_log)

Now get the trend/slope for each individual and make into a variable in df.

lizard_coeffs_CEWL <- data.frame(emtrends(lizard_CEWL_time_lm_log, 
                                     "lizard_ID", var = "time_min")) %>%
   left_join(lizard_tmt_info, 
             by = c("lizard_ID" = "individual_ID")) %>%
   dplyr::select(lizard_ID, rate_change = time_min.trend, 
                 SE, df, lower.CL, upper.CL,
                 date,
                 treatment,
                 mass_g, SVL_mm,
                 chamber_time_elapsed_hr
                 ) %>%
   dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs_CEWL)
##    lizard_ID   rate_change              SE                  df       
##  401    : 1   Min.   :-0.094984   Min.   :0.0005118   Min.   :24381  
##  402    : 1   1st Qu.:-0.033908   1st Qu.:0.0005524   1st Qu.:24381  
##  404    : 1   Median :-0.013033   Median :0.0005611   Median :24381  
##  406    : 1   Mean   :-0.014469   Mean   :0.0005875   Mean   :24381  
##  407    : 1   3rd Qu.: 0.004875   3rd Qu.:0.0005916   3rd Qu.:24381  
##  408    : 1   Max.   : 0.058490   Max.   :0.0009780   Max.   :24381  
##  (Other):26                                                          
##     lower.CL            upper.CL              date              treatment 
##  Min.   :-0.096089   Min.   :-0.093879   Min.   :2021-10-05   Control:11  
##  1st Qu.:-0.035009   1st Qu.:-0.032807   1st Qu.:2021-10-11   Cooling:11  
##  Median :-0.014197   Median :-0.011868   Median :2021-10-12   Heating:10  
##  Mean   :-0.015621   Mean   :-0.013318   Mean   :2021-10-13               
##  3rd Qu.: 0.003788   3rd Qu.: 0.005963   3rd Qu.:2021-10-18               
##  Max.   : 0.057419   Max.   : 0.059561   Max.   :2021-10-19               
##                                                                           
##      mass_g          SVL_mm      chamber_time_elapsed_hr     date_factor
##  Min.   : 8.80   Min.   :60.00   Min.   :0.9333          2021-10-05:6   
##  1st Qu.:10.38   1st Qu.:64.00   1st Qu.:1.0000          2021-10-11:6   
##  Median :11.85   Median :66.00   Median :1.0250          2021-10-12:7   
##  Mean   :11.68   Mean   :66.47   Mean   :1.1068          2021-10-18:7   
##  3rd Qu.:12.82   3rd Qu.:69.00   3rd Qu.:1.1792          2021-10-19:6   
##  Max.   :15.30   Max.   :75.00   Max.   :1.5333                         
## 
length(unique(lizard_coeffs_CEWL$lizard_ID))
## [1] 32

NOW we can look at how the RATE of CEWL change (CEWL change per MINUTE) is different among lizards.

change_LMM_1 <- lme4::lmer(data = lizard_coeffs_CEWL,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       mass_g + SVL_mm + 
                       chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))

drop1(change_LMM_1)
## boundary (singular) fit: see help('isSingular')
## Single term deletions
## 
## Model:
## rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr + 
##     (1 | date_factor)
##                         npar     AIC
## <none>                       -148.61
## treatment                  2 -120.58
## mass_g                     1 -150.42
## SVL_mm                     1 -150.33
## chamber_time_elapsed_hr    1 -150.41
anova(change_LMM_1)
## Analysis of Variance Table
##                         npar    Sum Sq   Mean Sq F value
## treatment                  2 0.0279063 0.0139532 39.5378
## mass_g                     1 0.0000263 0.0000263  0.0744
## SVL_mm                     1 0.0000813 0.0000813  0.2305
## chamber_time_elapsed_hr    1 0.0000538 0.0000538  0.1524

Drop mass first based on low SS and resulting AIC:

change_LMM_2 <- lme4::lmer(data = lizard_coeffs_CEWL,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       SVL_mm + 
                       chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))
drop1(change_LMM_2)
## boundary (singular) fit: see help('isSingular')
## Single term deletions
## 
## Model:
## rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr + 
##     (1 | date_factor)
##                         npar     AIC
## <none>                       -150.42
## treatment                  2 -121.58
## SVL_mm                     1 -152.32
## chamber_time_elapsed_hr    1 -152.19
anova(change_LMM_2)
## Analysis of Variance Table
##                         npar    Sum Sq   Mean Sq F value
## treatment                  2 0.0279308 0.0139654 39.8100
## SVL_mm                     1 0.0000111 0.0000111  0.0316
## chamber_time_elapsed_hr    1 0.0000680 0.0000680  0.1940

drop SVL:

change_LMM_3 <- lme4::lmer(data = lizard_coeffs_CEWL,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                             chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))
drop1(change_LMM_3)
## Single term deletions
## 
## Model:
## rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
##                         npar     AIC
## <none>                       -152.32
## treatment                  2 -123.20
## chamber_time_elapsed_hr    1 -154.15
anova(change_LMM_3)
## Analysis of Variance Table
##                         npar    Sum Sq   Mean Sq F value
## treatment                  2 0.0279241 0.0139620 41.4001
## chamber_time_elapsed_hr    1 0.0000513 0.0000513  0.1521

drop chamber_time_elapsed_hr:

change_LMM_4 <- lme4::lmer(data = lizard_coeffs_CEWL,
                          rate_change ~ 
                             treatment + 
                       (1|date_factor))
drop1(change_LMM_4)
## boundary (singular) fit: see help('isSingular')
## Single term deletions
## 
## Model:
## rate_change ~ treatment + (1 | date_factor)
##           npar     AIC
## <none>         -154.15
## treatment    2 -116.99
anova(change_LMM_4)
## Analysis of Variance Table
##           npar   Sum Sq  Mean Sq F value
## treatment    2 0.027917 0.013959  42.898

The simplest model is our best/final model.

Check assumptions:

plot(change_LMM_4)

Assumptions look good.

Save with p-values and remove the random effect because it did not account for any variation:

change_CEWL_best <- lm(data = lizard_coeffs_CEWL,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment)
summary(change_CEWL_best)
## 
## Call:
## lm(formula = rate_change ~ treatment, data = lizard_coeffs_CEWL)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.043019 -0.015382  0.001456  0.013347  0.037537 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.009175   0.006084  -1.508  0.14239    
## treatmentCooling -0.042790   0.008605  -4.973 2.73e-05 ***
## treatmentHeating  0.030128   0.008817   3.417  0.00189 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02018 on 29 degrees of freedom
## Multiple R-squared:  0.7057, Adjusted R-squared:  0.6854 
## F-statistic: 34.77 on 2 and 29 DF,  p-value: 1.981e-08

Treatment group explained a significant amount of the variation.

and get average change ~ tmt:

avg_CEWL_change_coeff <- data.frame(emmeans(change_CEWL_best, "treatment"))
avg_CEWL_change_coeff
##   treatment       emmean          SE df     lower.CL     upper.CL
## 1   Control -0.009174826 0.006084346 29 -0.021618711  0.003269059
## 2   Cooling -0.051965300 0.006084346 29 -0.064409185 -0.039521415
## 3   Heating  0.020952679 0.006381316 29  0.007901423  0.034003936
avg_CEWL_change_pairs <- data.frame(pairs(emmeans(change_CEWL_best, "treatment"))) # p values
avg_CEWL_change_pairs
##            contrast    estimate          SE df   t.ratio      p.value
## 1 Control - Cooling  0.04279047 0.008604565 29  4.972997 7.921430e-05
## 2 Control - Heating -0.03012750 0.008817055 29 -3.416958 5.223942e-03
## 3 Cooling - Heating -0.07291798 0.008817055 29 -8.270106 1.205048e-08

and finally, t-tests of rate of change per tmt group:

lizard_coeffs_CEWL_control <- lizard_coeffs_CEWL %>% 
  dplyr::filter(treatment == "Control")
t.test(lizard_coeffs_CEWL_control$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_CEWL_control$rate_change
## t = -2.64, df = 10, p-value = 0.02473
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.016918337 -0.001431314
## sample estimates:
##    mean of x 
## -0.009174826
lizard_coeffs_CEWL_cool <- lizard_coeffs_CEWL %>% 
  dplyr::filter(treatment == "Cooling")
t.test(lizard_coeffs_CEWL_cool$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_CEWL_cool$rate_change
## t = -6.9758, df = 10, p-value = 3.826e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.06856350 -0.03536711
## sample estimates:
##  mean of x 
## -0.0519653
lizard_coeffs_CEWL_heat <- lizard_coeffs_CEWL %>% 
  dplyr::filter(treatment == "Heating")
t.test(lizard_coeffs_CEWL_heat$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_CEWL_heat$rate_change
## t = 3.0047, df = 9, p-value = 0.01484
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.005178089 0.036727269
## sample estimates:
##  mean of x 
## 0.02095268

* CEWL ~ Temp

lizard_CEWL_temp_loglog_lm <- lm(data = temp_time_series, 
                                 log(CEWL_g_m2h) ~ log(calibrated_cloacal_temp) * lizard_ID)
plot(lizard_CEWL_temp_loglog_lm)

Now get the trend/slope for each individual and make into a variable in df.

lizard_coeffs_CEWL_temp <- data.frame(emtrends(lizard_CEWL_temp_loglog_lm, 
                                          "lizard_ID", var = "calibrated_cloacal_temp")) %>%
   left_join(lizard_tmt_info, 
             by = c("lizard_ID" = "individual_ID")) %>%
   dplyr::select(lizard_ID, rate_change = calibrated_cloacal_temp.trend, 
                 SE, df, lower.CL, upper.CL,
                 date,
                 treatment
                 ) %>%
   dplyr::mutate(date_factor = as.factor(date))
length(unique(lizard_coeffs_CEWL_temp$lizard_ID))
## [1] 32

then get average rate of change by treatment group:

avg_CEWL_temp_change <- lm(data = lizard_coeffs_CEWL_temp, rate_change ~ treatment)
summary(avg_CEWL_temp_change)
## 
## Call:
## lm(formula = rate_change ~ treatment, data = lizard_coeffs_CEWL_temp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42655 -0.01854 -0.00098  0.03850  0.33076 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.15748    0.04090  -3.850 0.000600 ***
## treatmentCooling  0.21314    0.05784   3.685 0.000935 ***
## treatmentHeating  0.17726    0.05927   2.991 0.005630 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1357 on 29 degrees of freedom
## Multiple R-squared:  0.3475, Adjusted R-squared:  0.3025 
## F-statistic: 7.722 on 2 and 29 DF,  p-value: 0.002048
avg_CEWL_temp_change_coeff <- data.frame(emmeans(avg_CEWL_temp_change, "treatment"))
avg_CEWL_temp_change_coeff
##   treatment      emmean         SE df    lower.CL    upper.CL
## 1   Control -0.15747877 0.04090155 29 -0.24113183 -0.07382571
## 2   Cooling  0.05566413 0.04090155 29 -0.02798893  0.13931719
## 3   Heating  0.01977732 0.04289790 29 -0.06795875  0.10751338
avg_CEWL_temp_change_pairs <- data.frame(pairs(emmeans(avg_CEWL_temp_change, "treatment")))
avg_CEWL_temp_change_pairs
##            contrast    estimate         SE df    t.ratio     p.value
## 1 Control - Cooling -0.21314290 0.05784352 29 -3.6848188 0.002611154
## 2 Control - Heating -0.17725609 0.05927197 29 -2.9905549 0.015087163
## 3 Cooling - Heating  0.03588681 0.05927197 29  0.6054601 0.818280857

and finally, t-tests of rate of change per tmt group:

lizard_coeffs_CEWL_temp_control <- lizard_coeffs_CEWL_temp %>% 
  dplyr::filter(treatment == "Control")
t.test(lizard_coeffs_CEWL_temp_control$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_CEWL_temp_control$rate_change
## t = -2.3009, df = 10, p-value = 0.04419
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.309978692 -0.004978853
## sample estimates:
##  mean of x 
## -0.1574788
lizard_coeffs_CEWL_temp_cool <- lizard_coeffs_CEWL_temp %>% 
  dplyr::filter(treatment == "Cooling")
t.test(lizard_coeffs_CEWL_temp_cool$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_CEWL_temp_cool$rate_change
## t = 4.905, df = 10, p-value = 0.0006185
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.03037812 0.08095015
## sample estimates:
##  mean of x 
## 0.05566413
lizard_coeffs_CEWL_temp_heat <- lizard_coeffs_CEWL_temp %>% 
  dplyr::filter(treatment == "Heating")
t.test(lizard_coeffs_CEWL_temp_heat$rate_change)
## 
##  One Sample t-test
## 
## data:  lizard_coeffs_CEWL_temp_heat$rate_change
## t = 2.8899, df = 9, p-value = 0.01789
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.004295961 0.035258674
## sample estimates:
##  mean of x 
## 0.01977732

Figures

COLOR

cool <- c(brewer.pal(11, "RdBu")[c(11)])
heat <- c(brewer.pal(11, "RdBu")[c(2)])
control <- c(brewer.pal(9, "Greys")[c(5)])
cntrls <- c(brewer.pal(9, "Greys")[c(7, 5, 3, 5, 4, 5, 3, 6, 7, 3, 6)])

cloacal temp ~ time

ggplot(temp_time_series) +
   aes(x = time_sec, 
       y = calibrated_cloacal_temp,
       group = lizard_ID,
       color = treatment) +
   geom_line(size = 0.2) +
  geom_smooth(se=F, method = "lm") +
   theme_classic() +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = "Internal Body Temperature (°C)", 
        color = "Treatment", se = FALSE) -> ctemp_time
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
ctemp_time
## `geom_smooth()` using formula = 'y ~ x'

ggplot(temp_time_series) +
   aes(x = log(time_sec), 
       y = log(calibrated_cloacal_temp),
       group = lizard_ID,
       color = treatment) +
   geom_line() +
   theme_classic() +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = "Internal Body Temperature (°C)", 
        color = "Treatment", se = FALSE) -> ctemp_time_log
ctemp_time_log

* cloacal temp (relative) ~ time

ggplot(temp_time_series) +
   aes(x = time_min, 
       y = relative_clo_temp,
       group = lizard_ID,
       color = treatment) +
  
   geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
  
   geom_line() +
   
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        plot.margin = margin(t = 1, r = 1, b = 1, l = 4.3, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_x_continuous(limits = c(1,15),
                     breaks = c(seq(1, 15, 2)),
                     labels = c(seq(1, 15, 2))) +
  scale_y_continuous(limits = c(-21,21),
                     breaks = c(-20, -10, 0, 10, 20),
                     labels = c(-20, -10, 0, 10, 20)) +
  
  annotate("text", x = 12, y = -2, label = "Control", parse = T,
           size = 2.9, family = "sans", color = control) +
  annotate("text", x = 4, y = -11, label = "Cooling", parse = T,
           size = 2.9, family = "sans", color = cool) +
  annotate("text", x = 4, y = 10, label = "Heating", parse = T,
           size = 2.9, family = "sans", color = heat) +
  
   labs(x = "Time (minutes)" , 
        y = expression('Relative T'[b]~' (°C)'),
        color = "Treatment", se = FALSE) -> ctemp_time_relative
ctemp_time_relative

* Temp Rate of Change

ggplot() +
  geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
  
  geom_point(data = lizard_coeffs_temp,
              aes(x = treatment,
                   y = rate_change,
                   shape = treatment,
                   color = treatment,
                   fill = treatment),
                   size = 1,
               alpha = 0.3,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = avg_temp_change_coeff,
                aes(x = treatment,
                    y = emmean, 
                    color = treatment,
                    group = treatment,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = 0.3,
                alpha = 0.9) +
  geom_point(data = avg_temp_change_coeff,
             aes(x = treatment,
                   y = emmean,
                 shape = treatment,
                 fill = treatment),
              color = "black",
              size = 3) +
  
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        plot.margin = margin(t = 1, r = 1, b = 1, l = 8, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  scale_fill_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
  ylim(-1.8, 1.8) +
  
   annotate("text", x = 1, y = 1.8, label = "a", size = 3) +
   annotate("text", x = 2, y = 1.8, label = "b", size = 3) +
   annotate("text", x = 3, y = 1.8, label = "c", size = 3) +
  
  annotate("text", x = 1, y = 1.55, label = "***", size = 4) +
  annotate("text", x = 2, y = 1.55, label = "***", size = 4) +
  annotate("text", x = 3, y = 1.55, label = "***", size = 4) +
  
   labs(x = "" , 
        y = expression(Delta ~ 'log-T'[b]~'  log-'*min^-1*''),
        color = "Treatment") -> rate_change_temp_fig
rate_change_temp_fig

# save
ggsave(filename = "temp_change_min.pdf",
       plot = rate_change_temp_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

dorsal temp ~ time

ggplot(temp_time_series) +
   aes(x = time_sec, 
       y = calibrated_dorsal_temp,
       group = lizard_ID,
       color = treatment) +
   geom_line() +
   theme_classic() +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x= "Time (seconds)" , 
        y= "Surface Body Temperature (°C)", 
        color ="Treatment") -> dtemp_time
dtemp_time

CEWL (raw) ~ time

ggplot(temp_time_series) +
   aes(x = time_sec, 
       y = CEWL_g_m2h,
       group = lizard_ID,
       color = treatment) +
   geom_line() +
   # geom_point(size = 0.01, 
   #            alpha = 0.25) +
   # geom_smooth(method = "loess", 
   #             size = 2.5, 
   #             se = F) +
   theme_classic() +
   theme(text = element_text(size = 15)) + 
   theme(axis.text.x = element_text(size = 10)) +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = bquote('CEWL (g/'*m^2*'h)'), 
        color = "Treatment", 
        se = FALSE) +
   scale_x_continuous(limits = c(0, 925), 
                      breaks=seq(0, 950, 100)) -> CEWL_time
CEWL_time

CEWL (MIN) ~ time

ggplot(by_minute_avg) +
   aes(x= time_min_block, 
       y = CEWL_min_avg,
       group = treatment,
       color = treatment
       ) +
   geom_point(size = 1, 
              alpha = 0.5) +
  geom_line() +
  geom_smooth(method = "lm",
              se = F) +
   theme_classic() +
   theme(text = element_text(size = 15)) + 
   theme(axis.text.x = element_text(size = 10)) +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = bquote('CEWL (g/'*m^2*'h)'), 
        color = "Treatment", 
        se = FALSE)  -> CEWL_time_min
CEWL_time_min
## `geom_smooth()` using formula = 'y ~ x'

legend

This code uses one of our plots to make a large, detailed legend, which we can save and use in ggarrange later.

ggplot(temp_time_series) +
   
   geom_smooth(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       color = treatment),
       method = "lm", 
        formula = y ~ poly(x, 2),
               size = 1, 
               se = F) +
  geom_point(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       fill = treatment,
       shape = treatment),
       size = 4, 
       alpha = 1) +
  
  scale_x_continuous(limits = c(14, 39),
                     breaks = c(seq(15, 35, 5)),
                     labels = c(seq(15, 35, 5))) +
  theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = c(0.5,0.5)) +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") +
   scale_fill_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment")  -> fancy_legend
fancy_legend # view

# save
LEGEND <- as_ggplot(get_legend(fancy_legend)) 

# export
ggsave(filename = "legend_only.pdf",
       plot = LEGEND,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 40, height = 30)

* CEWL (raw) ~ cloacal temperature

ggplot(temp_time_series) +
  aes(
    x = calibrated_cloacal_temp,
    y = CEWL_g_m2h,
    group = lizard_ID,
    color = treatment
  ) +
  geom_line() +

  # annotations for cooling group
   annotate(geom = "point", x = 13, y = 3, 
           size = 3, pch = 25, fill = cool) + 
   annotate("text", x = 11, y = 2.5, label = "T[15]", parse = T,
           size = 3, family = "sans", color = cool) +
   annotate(geom = "point", x = 36, y = 19, 
           size = 3, pch = 25, fill = cool) + 
   annotate("text", x = 38, y = 18.5, label = "T[1]", parse = T,
           size = 3, family = "sans", color = cool) +
  
  # annotations for heating group
   annotate(geom = "point", x = 15, y = 14, 
           size = 3, pch = 24, fill = heat) +  
   annotate("text", x = 13, y = 14, label = "T[1]", parse = T,
           size = 3, family = "sans", color = heat) +
   annotate(geom = "point", x = 38, y = 30, 
           size = 3, pch = 24, fill = heat) +
   annotate("text", x = 40, y = 30, label = "T[15]", parse = T,
           size = 3, family = "sans", color = heat) +
   
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
  scale_x_continuous(limits = c(10, 40),
                     breaks = c(seq(10, 40, 10)),
                     labels = c(seq(10, 40, 10))) +
  scale_y_continuous(limits = c(0, 31),
                     breaks = c(seq(0, 30, 10)),
                     labels = c(seq(0, 30, 10))) +
  
   labs(x = expression('T'[b]~' (°C)'),
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        se = FALSE) -> CEWL_ctemp
CEWL_ctemp

# save
ggsave(filename = "CEWL_int_body_temp.pdf",
       plot = CEWL_ctemp,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 150, height = 100)

CEWL (log-log) ~ cloacal temperature

ggplot(temp_time_series) +
  aes(x = log(calibrated_cloacal_temp), 
      y = log(CEWL_g_m2h),
      group = lizard_ID,
      shape = treatment,
      color = treatment) +
  geom_line() +
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "right",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   labs(x = "log- Internal Body Temperature (°C)" , 
        y = bquote('log- CEWL (g '*m^-2*' '*h^-1*')'), 
        se = FALSE) -> CEWL_ctemp_log
CEWL_ctemp_log

# save
ggsave(filename = "CEWL_int_body_temp.pdf",
       plot = CEWL_ctemp,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 150, height = 100)

CEWL (MIN) ~ cloacal temperature

ggplot(by_minute_indiv) +
  aes(
    x = log(clo_temp_per_min),
    y = log(CEWL_per_min),
    group = lizard_ID,
    shape = treatment,
    color = treatment
  ) +
  geom_point(size = 1, 
               alpha = 0.8) +
  geom_line() +
  # scale_x_continuous(limits = c(12.5, 40),
  #                    breaks = c(seq(15, 35, 5)),
  #                    labels = c(seq(15, 35, 5))) +

   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
   labs(x = "Internal Body Temperature (°C)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        se = FALSE)

CEWL (MIN avg) ~ cloacal temperature

ggplot(by_minute_avg) +
  aes(
    x = clo_temp_min_avg,
    y = CEWL_min_avg,
    #group = lizard_ID,
    shape = treatment,
    color = treatment
  ) +
   geom_point(size = 1, 
               alpha = 0.8) +
  geom_line() +
   geom_smooth(method = "lm",
        #formula = y ~ poly(x, 2),
               size = 1.5,
               se = F) +
  scale_x_continuous(limits = c(12.5, 40),
                     breaks = c(seq(15, 35, 5)),
                     labels = c(seq(15, 35, 5))) +

   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
   labs(x = "Internal Body Temperature (°C)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        se = FALSE)
## `geom_smooth()` using formula = 'y ~ x'

ggplot(by_temp_minute_avg) +
  aes(
    x = clo_temp_chunk,
    y = CEWL_min_avg,
    #group = lizard_ID,
    shape = treatment,
    color = treatment
  ) +
   geom_point(size = 1, 
               alpha = 0.8) +
  geom_line() +
   geom_smooth(method = "lm",
        #formula = y ~ poly(x, 2),
               size = 1.5,
               se = F)
## `geom_smooth()` using formula = 'y ~ x'

CEWL Control

# get control data only
temp_time_series %>%
  dplyr::filter(treatment == "Control") %>%
  # pipe into ggplot
  ggplot() +
   aes(x = calibrated_cloacal_temp, 
       y = CEWL_g_m2h,
       color = lizard_ID) +
  
   geom_point(size = 0.1, 
               alpha = 0.1) +
   geom_smooth(method = "lm", 
               size = 1, 
               se = F) +

   ylim(0,32) +
   xlim(25,31) +
  
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   scale_color_manual(values = cntrls) +
  
   labs(x = "Internal Body Temperature (°C)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        color = "Treatment", 
        se = FALSE) -> CEWL_control
CEWL_control
## `geom_smooth()` using formula = 'y ~ x'

# save
ggsave(filename = "CEWL_temp_control_lm.pdf",
       plot = CEWL_control,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)
## `geom_smooth()` using formula = 'y ~ x'

* CEWL (relative) ~ time

# check n lizard sample size
temp_time_series %>%
  dplyr::filter(complete.cases(time_min, relative_CEWL)) %>%
  group_by(lizard_ID) %>%
  summarise(max(relative_CEWL)) %>%
  nrow()
## [1] 32
# plot
ggplot(temp_time_series) +
  aes(
    x= time_min, 
    y = relative_CEWL,
    group = lizard_ID,
    color = treatment
  ) +
  
   geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
  
   geom_line() +
   
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        axis.title.y = element_text(vjust = 0.5),
        plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  
  scale_x_continuous(limits = c(1,15),
                     breaks = c(seq(1, 15, 2)),
                     labels = c(seq(1, 15, 2))) +
  scale_y_continuous(limits = c(-18,18),
                     breaks = c(-10, 0, 10),
                     labels = c(-10, 0, 10)) +
  
  # annotate("text", x = 12, y = -2, label = "Control", parse = T,      # too cluttered
  #          size = 5, family = "sans", color = control) +
  annotate("text", x = 4, y = -15, label = "Cooling", parse = T, #9.3
           size = 2.9, family = "sans", color = cool) +
  annotate("text", x = 4, y = 7, label = "Heating", parse = T,
           size = 2.9, family = "sans", color = heat) +
  
   labs(x = " ", #"Time (minutes)" , 
        y = bquote('Relative CEWL (g '*m^-2*' '*h^-1*')')) -> CEWLr_time
  
CEWLr_time

# save
ggsave(filename = "CEWL_relative_time.pdf",
       plot = CEWLr_time,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 150, height = 100)

* CEWL Rate of Change

ggplot() +
  geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
  
  geom_point(data = lizard_coeffs_CEWL,
              aes(x = treatment,
                   y = rate_change,
                   shape = treatment,
                   color = treatment,
                   fill = treatment),
                   size = 1,
               alpha = 0.3,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = avg_CEWL_change_coeff,
                aes(x = treatment,
                    y = emmean, 
                    color = treatment,
                    group = treatment,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .3,
                alpha = 0.9) +
  geom_point(data = avg_CEWL_change_coeff,
             aes(x = treatment,
                   y = emmean,
                 shape = treatment,
                 fill = treatment),
              color = "black",
              size = 3) +
  
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        axis.title.y = element_text(vjust = -1.8),
        plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  scale_fill_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
  scale_y_continuous(limits = c(-0.1, 0.1),
                     breaks = c(-0.1, 0, 0.1),
                     labels = c(-0.1, 0, 0.1)) +
  
   annotate("text", x = 1, y = 0.1, label = "a", size = 3) +
   annotate("text", x = 2, y = 0.1, label = "b", size = 3) +
   annotate("text", x = 3, y = 0.1, label = "c", size = 3) +
  
  annotate("text", x = 1, y = 0.086, label = "*", size = 4) +
  annotate("text", x = 2, y = 0.086, label = "***", size = 4) +
  annotate("text", x = 3, y = 0.086, label = "*", size = 4) +
  
   labs(x = "" , 
        y = expression(Delta ~ 'log-CEWL  log-'*min^-1*''),
        color = "Treatment") -> rate_change_CEWL_fig
rate_change_CEWL_fig

# save
ggsave(filename = "CEWL_change_min.pdf",
       plot = rate_change_CEWL_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

* CEWL-Temp Rate of Change

ggplot() +
  geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
  
  geom_point(data = lizard_coeffs_CEWL_temp,
              aes(x = treatment,
                   y = rate_change,
                   shape = treatment,
                   color = treatment,
                   fill = treatment),
                   size = 1,
               alpha = 0.3,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = avg_CEWL_temp_change_coeff,
                aes(x = treatment,
                    y = emmean, 
                    color = treatment,
                    group = treatment,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .3,
                alpha = 0.9) +
  geom_point(data = avg_CEWL_temp_change_coeff,
             aes(x = treatment,
                   y = emmean,
                 shape = treatment,
                 fill = treatment),
              color = "black",
              size = 3) +
  
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        axis.title.y = element_text(vjust = -0.5),
        plot.margin = margin(t = 1, r = 1, b = 4, l = 1, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  scale_fill_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
  scale_y_continuous(limits = c(-0.6, 0.6),
                     breaks = c(-0.5, 0, 0.5),
                     labels = c(-0.5, 0, 0.5)) +
  
  # pairwise comparisons between tmt groups
   annotate("text", x = 1, y = 0.6, label = "a", size = 3) +
   annotate("text", x = 2, y = 0.6, label = "b", size = 3) +
   annotate("text", x = 3, y = 0.6, label = "b", size = 3) +
  
  # t-tests for each tmt group
  annotate("text", x = 1, y = 0.52, label = "*", size = 4) +
  annotate("text", x = 2, y = 0.52, label = "**", size = 4) +
  annotate("text", x = 3, y = 0.52, label = "*", size = 4) +
  
   labs(x = "" ,
        y = expression(Delta ~ 'log-CEWL  log-T'[b]),
        color = "Treatment") -> rate_change_CEWL_temp_fig
rate_change_CEWL_temp_fig

# save
ggsave(filename = "CEWL_temp_change.pdf",
       plot = rate_change_CEWL_temp_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

* Experimental Design Skematic (fig 1)

We created reference values to use to make a figure that would show each step of our experiment.

exp_design <- read.csv("./Data/exp design.csv", 
                            header = TRUE,
                            fileEncoding="UTF-8-BOM")
ggplot(exp_design) +
   aes(x = time_min, 
       y = temp_C,
       color = treatment) +
     geom_line(size = 1) +
   geom_vline(xintercept = 60,
              #alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(30, 38, family = "sans"), label = "incubation (a)", 
             color = "black", size = 2.5) +
   geom_vline(xintercept = 65,
              #alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(62.5, 38, family = "sans"), label = "(b)", 
             color = "black", size = 2.5) +
   geom_vline(xintercept = 80,
              #alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(72.5, 38, family = "sans"), label = "tmt (c)", 
             color = "black", size = 2.5) +
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        axis.text.x = element_blank(),
        #axis.line.x = element_blank(),
        axis.ticks.x = element_blank()) +

   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time" , 
        y = expression('Target T'[b]~' (°C)'),
        color = "Treatment") +
   
    annotate(geom = "point", x = 80, y = 25, 
           size = 3, pch = 21, fill = control) +
   annotate(geom = "point", x = 80, y = 20, 
           size = 3, pch = 25, fill = cool) + 
   annotate(geom = "point", x = 80, y = 30, 
           size = 3, pch = 24, fill = heat) +
  
  annotate("text", x = 35, y = 24, label = "Control", parse = T,
           size = 2.9, family = "sans", color = control) +
  annotate("text", x = 18, y = 30, label = "Cooling", parse = T,
           size = 2.9, family = "sans", color = cool) +
  annotate("text", x = 18, y = 20, label = "Heating", parse = T,
           size = 2.9, family = "sans", color = heat) +
   
   scale_y_continuous(limits = c(15, 38), 
                      breaks=seq(15, 35, 5)) -> methods_schmatic_fig
methods_schmatic_fig

# save
ggsave(filename = "methods_schematic.pdf",
       plot = methods_schmatic_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

Save Grouped Figures

* CEWL & Tb x relative & rate of change (Fig 2)

CEWL_temp_multi <- ggarrange(
  CEWLr_time, rate_change_CEWL_fig, ctemp_time_relative, rate_change_temp_fig,
  nrow = 2, ncol = 2,
  widths = c(2,1.2),
  labels = c("A", "B", "C", "D"),
  label.x = 0.01,
  label.y = 1.015
)

ggsave(filename = "CEWL_temp_relative_and_change.pdf",
       plot = CEWL_temp_multi,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 160, height = 140)

* CEWL ~ Tb

CEWL_temp_relation <- ggarrange(
  CEWL_ctemp, rate_change_CEWL_temp_fig,
  nrow = 1, ncol = 2,
  widths = c(2,1.2),
  labels = c("A", "B"),
  #label.x = 0.01,
  label.y = 1.015
)

ggsave(filename = "CEWL_versus_temp.pdf",
       plot = CEWL_temp_relation,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 160, height = 70)

Relative CEWL & Rate of Change

multi_fig_relative <- ggarrange(CEWLr_time, 
          ggarrange(rate_change_CEWL_fig, LEGEND, 
                    nrow = 1, ncol = 2,
                    widths = c(2,1),
                    align = "hv"
                    ),
          nrow = 2, ncol = 1,
          heights = c(2,1),
          labels = c("A", "B"))
multi_fig_relative

ggsave(filename = "CEWL_time_multi_fig.pdf",
       plot = multi_fig_relative,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 120, height = 160)

Raw CEWL ~ Cloacal temp & control zoom-in

multi_fig_CEWL_by_temp <- ggarrange(CEWL_ctemp, 
                  ggarrange(CEWL_control, LEGEND, 
                            nrow = 1, ncol = 2,
                            widths = c(2, 1),
                            align = "hv"
                            ),
                  nrow = 2, ncol = 1,
                  heights = c(2,1),
          labels = c("A", "B"))
## `geom_smooth()` using formula = 'y ~ x'
multi_fig_CEWL_by_temp

ggsave(filename = "CEWL_temp_multi_fig.pdf",
       plot = multi_fig_CEWL_by_temp,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 120, height = 160)